Monte Carlo results for the hydrogen Hugoniot 
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Abstract 

We propose a theoretical Hugoniot obtained by combining results for the equation of state (EOS) 
from the Direct Path Integral Monte Carlo technique (DPIMC) and those from Reaction Ensemble 
Monte Carlo (REMC) simulations. The main idea of such proposal is based on the fact that DPMIC 
provides first-principle results for a wide range of densities and temperatures including the region 
of partially ionized plasmas. On the other hand, for lower temperatures where the formation of 
molecules becomes dominant, DPIMC simulations become cumbersome and inefficient. For this 
region it is possible to use accurate REMC simulations where bound states (molecules) are treated 
on the Born-Oppenheimer level using a binding potential calculated by Kolos and Wolniewicz. The 
remaining interaction is then reduced to the scattering between neutral particles which is reliably 
treated classically applying effective potentials. The resulting Hugoniot is located between the 
experimental values of Knudson et al. Q] and Collins et al. Q|. 

PACS numbers: 64.30,+t, 05.30.-d, 62.50.+p 
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The H-plasma is a very important and interesting many particle system. Hydrogen is the 
simplest and at the same time the most abundant element in the universe. Due to its high 
relevance for modern astrophysics, inertial confinement fusion and fundamental understand- 



ing of condensed matter, hydrogen continues to 
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y studied both, experimentally 



141 ]. At high temperatures and 



pressures, the hydrogen behavior is defined by the interaction between free electrons and 
protons (plasma state). With decreasing temperature, the contribution of bound states such 
as atoms and molecules to the EOS of hydrogen becomes of increasing importance, which 
at low temperatures completely define the hydrogen properties. Furthermore, as pointed 
out in many papers (Norman and Starostin ^> Ebeling et al. [la ]. Haronska et al. jl7|. 
Saumon and Chabrier 18]) there are strong theoretical arguments for a phase transition 
between two plasma phases. This issue which is of importance, for example, for models of 
Jovian planets is still actively debated. Among other important questions we mention the 
high-pressure compressibility, details of the pressure ionization and dissociation. 

For this reason, in the last decades considerable experimental and theoretical investi- 
gations were carried out to accurately determine the EOS of hydrogen at high pressures. 
Experimentally, the EOS for this region can be obtained using shock-wave techniques. The 
results of these experiments are usually discussed in form of an Hugoniot 



E = E + hp + Po ) (--- 
2 \p Po, 



(1) 



where the specific internal energy E at a state with the density p and the pressure p is 
connected to the initial conditions with the density po, the pressure po and the internal 
energy E . 

One of the well established experimental techniques for the creation of shock waves uses 
gas gun devices. With gas gun experiments, Nellis et al. reached maximum pressures of 20 
GPa and temperatures of 7000 K. More advanced techniques, the laser-driven experiments 
used by Collins et al. P and Da Silva et al. Q, allow to reach pressures up to 300 
GPa. At such pressures, as expected, hydrogen transforms from a molecular to a metallic 
state jf| • The results of laser-driven experiments have shown an unusual high compression 
p/po = 6 of deuterium, which deviates si gnifi cantly from a maximum compression of p/po = 
4.4 obtained within the SESAME EOS However, the experiments of Knudson et al. 

fi 

[l| which used magnetically driven flyer techniques (Z-pinch) do not support such high 
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TABLE I: Thermodynamic properties of deuterium plasma calculated by DPIMC 
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compressibilities and are close to those of SESAME and Restricted Path Integral Monte 
Carlo (RPIMC) []| results. The reason for this discrepancy of the two experiments is not 
yet completely understood and requires more detailed study j^j], including independent 
theoretical investigations which is the aim of this paper. It is also necessary to mention 
other important experimental techniques such as the convergent geometry technique 
The experimental point obtained by Belov et al. within this technique is located between 
the results of laser — driven and magnetically driven flyer experiments. 

An Hugoniot can be also determined theoretically from the equation of state. This enables 
us to compare different theoretical approaches and computer simulations with experimental 
results, which cover a large region in the phase diagram of hydrogen. They start at tem- 
peratures of about 20 K and at a density of p = 0.171 g/cm 3 , which corresponds to the 
liquid state, and go up to temperatures and densities where only free electrons and nuclei 
exist. To our knowledge, there is no theory or computer simulation which rigorously and 
consistently describes the complete region of the EOS achievable by experiments. For ex- 
ample, the linear mixing model (LM) of Ross 2l( rather well predicts the behavior of the 
laser driven experiments; however it is a semi-empirical theory which interpolates between 
molecular and metallic states of hydrogen. 



Further, the region of completely and partially ionized hydrogen can be described ana 
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23j . In such methods, an EOS is 
obtained from a fugacity expansion (ACTEX) 23j and modified fugacity expansions which 



lytically using the methods of quantum statistics 
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are upgraded by means of quantum-field theoretical methods (leading to dynamical screen- 



ing, self energy and lowering of the ionization energy 
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to condense the results in form of Pade approximations 



221). In the latter case it is useful 



24( | . (from Debye to Gellman- 



Brueckner). Of course, the EOS following from these theories cannot reproduce the Hugo- 
niot over the full range of density and pressure. It gives only the asymptotic behavior at 
higher temperatures. The typical behavior of the analytical theory [24] is shown in Fig. 
It coincides only asymptotically with the ab initio RPIMC calculations and, with decreasing 
temperature, deviates considerably from those results. The Hugoniot calculated within the 
ACTEX theory which is not shown here exhibits a similar behavior 2^ |. 

The main reason for the failure of the analytical theories is obvious. As we mentioned al- 
ready, for lower temperatures, the neutral particles, i.e., H-atoms and H 2 — molecules, become 
increasingly important, giving rise to a strongly coupled dense gas or liquid. Under such 
conditions it is necessary to invoke the methods of the theory of liquids. The simplest th eory 
for this region is the fugacity expansion of the EOS up to the second virial coefficient [25]. 
This theory is applicable only for low densities and cannot correctly describe the molecular 
dissociation which is an important process occurring during shock wave experiments. For 
moderate densities, the fluid variational theory (FVT), proposed by Ross et al. 26], can be 
applied. This theory was further developed by Juranek and Redmer Q] to many component 
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FIG. 1: Experimental and theoretical results for the deuterium Hugoniot 
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systems, where molecular dissociation occurs. The effective interaction potentials |26l. l27| 
between components used within FVT achieve good agreement with experimental gas gun 
data of Nellis et al. |3], Fig. □ 

A powerful tool for the investigation of the hydrogen EOS is ab initio computer simulation. 
Quantum molecular dynamics simulations, based on a density functional theory, are usually 
applied to investigate the atomic and molecular region j^, but have difficulties to describe 
the partially ionized plasma. The wave packet molecular dynamics also covers the region of 
the fully ionized plasma but yields unexpectedly high compressibilties. In this work we 
will not discuss these methods in detail and refer to the work cited. 

The Path Integral Monte Carlo method is another first principle method which is well 
suited for the investigation of the EOS of hydrogen over a wide parameter range. Except 
for the problem of the Fermi statistics, it is an exact solution of the many-body quantum 
problem for a finite system in thermodynamic equilibrium. The treatment of the "sign 
problem" makes the main difference between the RPIMC method used by Ceperley and 
Militzer [3 and the Direct Path Integral Monte Carlo (DPIMC) method used by Filinov et 
al. Q, 13 an d others. This problem is beyond the present paper, here we restrict ourselves 
to the discussion of the DPIMC method. 

The idea of DPIMC is the well known: any thermodynamic property of a two-component 
plasma with N e electrons and N p protons at a temperature T and volume V is defined by 
the partition function Z(N e , N p , V, T): 

Z{N ei N P1 V, T) = Y, j y d 1 dr Pfa r > °\ n (2) 

where q (r) comprises the coordinates of the protons (electrons), a stands for the spin of 
the electrons, and p is the density matrix of the system. Taking into account the electron 
spin and the Fermi statistics (antisymmetrization), the density matrix is expressed by a 
path integral where all electrons are represented by fermionic loops with a number of 
intermediate coordinates (beads). In our simulations, we used an effective quantum pair po- 
tential, which is finite at zero distance 29] . This potential was obtained by Kelbg as a result 
of a first-order perturbation theory. The simulation has been performed at temperatures of 
10 4 K and higher in a wide range of particle densities. Under these conditions the exchange 
effects for protons are negligible. In the present calculations, we used an improved treatment 
of the electron exchange, i.e., we took into account the exchange interaction of electrons from 



TABLE II: Hugoniot data calculated by REMC 
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neighbor Monte Carlo cells, namely first from the nearest neighbor cells (3 3 — 1), then from 
the next neighbors (5 3 — 1) and so on. The calculated thermodynamic properties of hydrogen 
allowed us to compute the shock Hugoniot of deuterium using Eq. (jHJ) 

H = E-E - ±(p+ Po )(V - V ) = 0. (3) 

Following the work Q| we chose p = 0, p — 0.171 g/cm 3 , E = —15.886 eV per atom and 
computed the pressure pi and the energy Ei at a given constant temperature T (from 10 K 
to 10 6 K) and three values of the volume Vi = l/pi corresponding to r s =1.7, 1.86, and 2, 
where r s = f/as, f = (3/47m p ) 1//3 , n p is the particle density, a B - Bohr radius. The results of 
the calculations are presented in Table 1. Substituting the obtained values pi, Ei and Vi into 
the Hugoniot we determine the volume range V\, V-i where the function H(p, V, E) changes 
its sign. The value of the density at the Hugoniot is calculated by linear interpolation of 
the function H between V\ and Vi- The values of the pressure and of the total energy are 
shown in the Table 1 only for those density values between which the value of the density 
lies on the Hugoniot at a given temperature. The values of density and pressure on the 
Hugoniot are placed in the last two columns of Table 1. and are polotted together with 
selected theoretical and experimental data in Fig. ^ The lowest temperature included in 
this figure for the DPIMC is 15625 K. 

In order to correctly describe the quantum mechanics of the formation of molecules at 
temperatures lower than 10000 K, it is necessary to take many beads. In this region, DPIMC 
calculations become very time consuming and the convergence is poor. The natural proposal 
which appears for this region is to use the asymptotic property of the path integral which, 
for heavy particles, goes over into the classical partition function. For such systems, the 
classical Monte Carlo scheme can be applied. An advanced version of the classical Monte 
Carlo scheme is the reaction ensemble Monte Carlo technique (REMC) [3j|. This method 
incorporates the quantum mechanical description of bound states, while the scattering states 
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FIG. 2: Results for the combined Hugoniot 

are treated classically. As was shown by Bezkrovniy et al. [14], REMC describes the low 
temperature region very well, and yields good agreement with the gas gun experiments by 
Nellis et al. |3( Fig. ^ In these simulations the energy levels for the molecular partition 
functions of hydrogen and deuterium are obtained by solving the Schrodinger equation with 
the potential calculated by Kolos and Wolniewicz |31[. On the basis of the REMC, results 
are obtained much easier as compared to those from molecular dynamics based on density 
functional theory; see Bonev et al. [l3[ and Fig. 2. Our REMC data are presented in Table 
2. 

In order to get a unified picture combining DPIMC and REMC, we use the fact that 
REMC turns out to be the limiting case of DPIMC at low temperatures, where hydrogen 
consists only of atoms and molecules. Therefore, it is obvious to use the asymptotic results 
of both methods to construct an Hugoniot which can be applied in the entire range of 
compression. For the construction of the combined Hugoniot we carefully analyzed the 
region where the Hugoniots produced by the two methods can be connected to each other. 
As we can see from Fig. Qthe Hugoniot calculated within DPIMC ends at the point 15625 K. 
At this temperature, the largest contribution to the EOS are given by molecular states. As 
natural continuation of the DPIMC Hugoniot, we take the point of 15000 K produced by 
REMC. We want to stress here that these two methods are completely independent and 
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no interpolation procedure is used. Just two points at 15625 K of DPIMC and 15000 K 
of REMC are connected to each other. The final Hugoniot is plotted in Fig. |21 and shows 
a maximum compressibility of approximately 4.75 as compared to the initial deuterium 
density. 
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